function sim = mit_shock_entrants(sim, eq_SS, param, glob, options)

% Solves "backwards" and iterates "forwards" on one iteration of the MIT
% shock

% start at the end and solve backwards
v2    = reshape(eq_SS.v.v2, glob.n(1), glob.n(2)); % assume that the expected value is in SS at date T

tic() 

for t = options.T:-1:1
    
    param.C = sim.C(t);
    param.A = sim.A(t); 
    param.w = sim.W(t);
    
    
    param.mu_f = sim.cF(t);
    param.ce = sim.cE(t);
    
 %   v2    = griddedInterpolant(glob.B, glob.S,v2, 'linear');   
    
    %% find D that sets value of entry equal to 0
    tosolve = @(D) solve_b_minimal(D, v2, param, glob, options);
   
    if t == 1
        disp('t = 1')
    end
    
    %D = fzero(tosolve, eq_SS.D);
    try
        D = bisect(tosolve, .5*eq_SS.D, 2*eq_SS.D);
    catch
        disp('D not bracketed')
    end
        
    param.D = D;
    
    [~,ve, v] = solve_b(D, v2, param, glob, options);
    sim.ve(t) = ve.eve;
    sim.D(t) = D;
    % update v2
    v2    = reshape(v.v2, glob.n(1), glob.n(2));

    %% continue with incumbent problem
%     % solve on finer grid
%     tmp      = glob.PnK;
%     glob.PnK = glob.PnKf;
%     vf       = solve_valfunc([v2; v2],glob.sf,param,glob,options);
%     glob.PnK = tmp;
    
    vf = v;
    
    % Record optimal policy
    sim.y(:,t)  = vf.y; 
    sim.l(:,t)  = vf.l;
    sim.p(:,t)  = vf.p;  
    
    sim.exit(:,t) = vf.p_exit;
        
end
toc()

%% Now iterate forwards on policies
tic()

fspaceerg   = fundef({'spli',glob.bgridf,0,1});


for t = 1:options.T
    
    % Distribution from last period
   if t > 1
       L               = sim.mu(:,t-1);
       
       lp = sim.l(:,t-1);
       le = glob.ve_l;
       
       p_exit          = sim.exit(:,t-1);
   else
        L               = eq_SS.L;
        
        lp = eq_SS.l;
        le = eq_SS.ve_l;
        
        p_exit          = eq_SS.p_exit;
   end
       
    % distribution of entrants
    G                  = zeros(glob.Nsf,1);
    G((glob.sf(:,1) == glob.bgridf(end))) = glob.a_stat_E;

    entrant_policy = repmat(le, 1, glob.nf(1))';
    entrant_policy = entrant_policy(:)';
    
    entrant_policy          = min(entrant_policy, glob.maxb);
    entrant_policy          = max(entrant_policy, glob.minb);

    Qb             = funbas(fspaceerg, entrant_policy');

    Qent           = dprod(glob.QeI,Qb);
    G = Qent'*G;
   
    % set up transition matrix
    lp          = min(lp, glob.maxb);
    lp          = max(lp, glob.minb);
      
    Qb              = funbas(fspaceerg, lp);
    Q               = dprod(glob.Qef,Qb);
        
     param.D = sim.D(t);
     param.C = sim.C(t);
    param.w = sim.W(t);
    yf = sim.y(:,t);
    
    tosolve  = @(M) iterate_dist(M, L, G, Q, p_exit,yf, param,glob,options);
    
    try
        M        = bisect(tosolve, max(0,eq_SS.M - 1), eq_SS.M + 5);
    catch
        M = 0; 
    end
        
    sim.M(t) = M;
    
    
    sim.mu(:,t) = Q'*((1-param.gamma).*(1-p_exit).*L) + M*G;

    sim.entry_mass(t) = M;
    sim.entry_rate(t) = sim.entry_mass(t)/sum(sim.mu(:,t));

    eq_mit.L          = sim.mu(:,t);
    eq_mit.yf         = sim.y(:,t);
    
    aggs              = find_agg_mit(eq_mit, param, glob, options);
    
    sim.C(t)          = aggs.C;        
                                  
    L                 = sum(eq_mit.L.*sim.l(:,t)); % this is actually the previous period's labor supply - FIXED WLG 7/30                    
    sim.L(t)          = L;
    sim.W(t)          = L^param.nu*param.psi;
      
    % get entrant info
        % distribution of firms ages a=0,1,2,3,4,5
    sim.a0(:,t)      = M*G;
    
    if t == 1
        sim.a1(:,t)      = Q'*((1-param.gamma)*(1-p_exit).*sim.a0(:,t));
        sim.a2(:,t)      = Q'*((1-param.gamma)*(1-p_exit).*sim.a1(:,t));
        sim.a3(:,t)      = Q'*((1-param.gamma)*(1-p_exit).*sim.a2(:,t));
        sim.a4(:,t)      = Q'*((1-param.gamma)*(1-p_exit).*sim.a3(:,t));
        sim.a5(:,t)      = Q'*((1-param.gamma)*(1-p_exit).*sim.a4(:,t));
    else  
        sim.a1(:,t)      = Q'*((1-param.gamma)*(1-sim.exit(:,t-1)).*sim.a0(:, t-1));
        sim.a2(:,t)      = Q'*((1-param.gamma)*(1-sim.exit(:,t-1)).*sim.a1(:, t-1));
        sim.a3(:,t)      = Q'*((1-param.gamma)*(1-sim.exit(:,t-1)).*sim.a2(:, t-1));
        sim.a4(:,t)      = Q'*((1-param.gamma)*(1-sim.exit(:,t-1)).*sim.a3(:, t-1));
        sim.a5(:,t)      = Q'*((1-param.gamma)*(1-sim.exit(:,t-1)).*sim.a4(:, t-1));
    end
    
    if t == 6
        disp('hello, it is t = 6')
    end
    
    sim.employment_a0(t) = sum(sim.l(:,t).*sim.a0(:,t));
    sim.employment_a1(t) = sum(sim.l(:,t).*sim.a1(:,t));
    sim.employment_a2(t) = sum(sim.l(:,t).*sim.a2(:,t));
    sim.employment_a3(t) = sum(sim.l(:,t).*sim.a3(:,t));
    sim.employment_a4(t) = sum(sim.l(:,t).*sim.a4(:,t));
    sim.employment_a5(t) = sum(sim.l(:,t).*sim.a5(:,t));

    sim.employment_young(t) = sim.employment_a0(t) + sim.employment_a1(t) + ...
                              sim.employment_a2(t) + sim.employment_a3(t) + ...
                              sim.employment_a4(t) + sim.employment_a5(t);
    
                              % employment at firms of age < 1
    sim.entrant_employment(t) = sim.employment_a0(t);

                          
    rev = sim.p(:,t).*sim.y(:,t);
                          
    sim.output_a0(t) = sum(rev.*sim.a0(:,t));
    sim.output_a1(t) = sum(rev.*sim.a1(:,t));
    sim.output_a2(t) = sum(rev.*sim.a2(:,t));
    sim.output_a3(t) = sum(rev.*sim.a3(:,t));
    sim.output_a4(t) = sum(rev.*sim.a4(:,t));
    sim.output_a5(t) = sum(rev.*sim.a5(:,t));
    
    sim.output_young(t) = sim.output_a0(t) + sim.output_a1(t) + ...
                              sim.output_a2(t) + sim.output_a3(t) + ...
                              sim.output_a4(t) + sim.output_a5(t);
    
end
toc()



end